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ABSTRACT 

Aims. We compare stellar catalogues with position and proper motion components using a decomposition on a set of orthogonal vector 
spherical harmonics. We aim to show the theoretical and practical advantages of this technique as a result of invariance properties and 
the independence of the decomposition from a prior model. 

Methods. We describe the mathematical principles used to perform the spectral decomposition, evaluate the level of significance of 
the multipolar components, and examine the transformation properties under space rotation. 

Results. The principles are illustrated with a characterisation of systematic effects in the FK5 catalogue compared to Hipparcos and 
with an application to extraction of the rotation and dipole acceleration in the astrometric solution of QSOs expected from Gaia. 
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1. Introduction 

The differences between the positions of a set of common sources in two astrometric catalogues are conveniently described math- 
ematically by a vector field on a sphere. Each vector materialises the difference between the two unit vectors giving the direction 
of the common sources in each catalogue. This feature extends to the differences between the proper motions of each source which 
also generates a spherical vector field. Typically when connecting two positional catalogues with n common sources, one has the 
coordinates a' 1 ',^- 1 ' for the z'th source in the first catalogue and ctf* ,5^ for the same source in the second catalogue. Provided the 
two catalogues are close to each other, the differences can be mapped as a vector field with components in the local tangent plane 
given by 

V = [Aacos<5 = (ar (2) -ar (1) )coS(5 (1) , A6 = 5 {2) - 5 m ] (1) 

for each common source between the two catalogues. We wish to analyse these fields in order to summarise their largest overall or 
local features by means of a small set of base functions, which is much smaller in any case than the number of sources. 

Several papers in the astronomical literature have applied this overall idea with either scalar or vectorial functions. As far as 
we have been able to trace the relevant publications regarding the use of vector spherical harmonics in this context, this has been 
initiated by Mignard & Morando ( 1990 ) and published in a proceedings paper that is not widely accessible. One of the motivations 
of this paper is therefore to update and expand on this earlier publication and to provide more technical details on the methodology. 

The general idea of the decomposition has its root in the classical expansion of a scalar function defined on the unit sphere with 
a set of mutually orthogonal spherical harmonics functions. The most common case in astronomy and geodesy is the expansion 
of the gravitational potential of celestial bodies, in particular that of the Earth. The generalisation of such expansions to vectors, 
tensors, or even arbitrary fields was introduced in mathematics and theoretical physics decades ago (Gelfand, Minlos & Shapiro 
1963). For example, these generalisations are widely used in gravitational physics ( Throne| 1 980| |Suen| 1 986| |Blanchet & Damour| 



1986). Restricting it to the case of vectors fields, which is the primary purpose of this paper, we look for a set of base functions 
that would allow representing any vector field on the unit sphere as an infinite sum of fields, so that the angular resolution would 
increase with the degree of the representation (smaller details being described by the higher harmonics). One would obviously like 
for the lowest degrees to represent the most conspicuous features seen in the field, such as a rotation about any axis or systematic 
distortion on a large scale. It is important to stress at this point that expanding a vector field on this basis function is not at all the 
same as expanding the two scalar fields formed by the components of the vector field: the latter would depend very much on the 
coordinate system used, while a direct representation on a vectorial basis function reveals the true geometric properties of the field, 
regardless of the coordinate system, in the same way as the vector field can be plotted on the sphere without reference to a particular 
frame without a graticule drawn in the background. 
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Several related methods using spherical functions to model the errors in astrometric catalogues, the angular distribution of proper 
motions, or more simply as a means of isolating a rotation have been published and sometimes with powerful algorithms. |Brosche] 



(1966 1970 1 was probably the first to suggest the use of orthogonal functions to characterise the errors in all-sky astrometric 



catalogue, and his method was later improved by |Schwan| ([1977} to allow for a magnitude dependence. But in both cases the 
analysis was carried out separately for each component Aa cos 5, A6 of the error, by expanding two scalar functions. Therefore, this 
was not strictly an analysis of a vector field, but that of its components on a particular reference frame. The results therefore did 
not describe the intrinsic geometric properties of the field, which should reveal properties not connected to a particular frame (see 
also Appendix |Cj. Similarly, a powerful method for exhibiting primarily the rotation has been developed in a series of papers of 



Vityazev (2010 and references therein). The vector spherical harmonics (hereafter abridged in VSH) have been in particular used 



in several analyses of systematic effects in VLBI catalogues or in comparison of reference frames, such as in Arias et al. (2000), 
Titov & Malkin (20 09), or Gwinn et al. ( 1997[ l, the galactic velocity field (Makarov & Murp hy~[2007[ l, the analysis of zonal errors in 



space astrometry (Makar ov et al.|2012[ ), but none of these papers deals with the fundamental principles, the very valuable invariance 
properties and the relevant numerical methods. Within the Gaia preparatory activities the VSH are used also to compare sphere 
solutions and in the search of systematic differences at different scales by Bucciarelli et al ( 201 l| l and more is planned in the future 
to characterise and evaluate the global astrometric solution. 

In this paper we provide the necessary theoretical background to introduce the VSHs and the practical formulas needed to 
explicitly compute the expansion of a vector field. The mathematical principles are given in Section [2] with a few illustrations to 
show the harmonics of lower degrees. In the next section, Section [3] we emphasise the transformation of the VSHs and that of the 
expansion of a vector field under a rotation of the reference frame, with applications to the most common astronomical frames. 
Then in Section|4]we discuss the physical interpretation of the harmonics of first degree as the way of representing the global effect 
shown by a vector field, such as the axial rotation and its orthogonal counterpart, for which we have coined the term glide, and show 
its relationship to the dipole axial acceleration resulting from the Galactic aberration. The practical implementation is taken up in 
Section [5] where we also discuss the statistical testing of the power found in each harmonic. The results of particular applications 
to the FK5 Catalogue and to the simulated QSO catalogue expected from Gaia are respectively given in Sects. [6] and [7] Appendix 
|A]provides the explicit expressions of the VSH up to degree I = 4, while Appendix [B] deals with the practical formulas needed for 
numerical evaluation of the VSH, and Appendix [C] focuses on the relationship between expansions of the components of a vector 
field on the scalar spherical harmonics and the expansion of the same field on the VSHs. 



2. Mathematical principles 




Fig. 1. Local frame associated to the spherical coordinates a, 6, with the unit vectors along the longitude and latitude lines. 
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2.1. Definitions 



In this section we give the main definitions needed to use the VSHs and the theoretical background from which practical numerical 
methods are established and discussed later in the paper. The VSHs form an orthogonal set of basis functions for a vector field on a 
sphere and split into two categories referred to as the toroidal T\ m and spheroidal S/ m functions (in the physical literature the former 
are also called "magnetic" or "stream", and the latter "poloidal", "potential", or "electric"): 



1 



1) 

Tim — —U X Sl m . 



rVY lm = uxT 



Im • 



(2) 
(3) 



where u = r/r, r = \r\, r is the radius-vector of the point, and V denotes the gradient operator. Clearly, for points on the surface of 
a unit sphere r — 1 . It is also obvious that at each point on the sphere and for each I and m one has u ■ Si m = 0, u ■ T/ m = 0, and 
Tim ■ S[ m = 0. Taking into account that (see Fig. [TJ 



( cos a cos 5 \ 
sin a cos 5 
sin 6 
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cos 6 da 
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■ cos of sin 6 ^ 
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(so that \u\ = 1, \e a \ = 1, \e$\ = 1), one gets explicit formulas 



Ti m {a, 8) = 



1 



V7(Z + 1) 



85 



1 dY, 



im 



cos 5 da 



for the toroidal functions, and 



S, m (a, 5) 



1 



+ 1) 



1 dY lm dY lm 
-e„ + — —e s 



cos 5 da 



86 



for the spheroidal functions. The Y\ m are the standard spherical functions defined here with the following sign convention 



YUa, 5) = (-!)" 



21+1 (Z-m)! 



P /m (sin<5) e " 



, An (1 + m)\ 
for m > and 

Y l - m (a,5) = (-lTY; m (a,5) 

for m < 0, where superscript '*' denotes complex conjugation. The associated Legendre functions are defined as 



PUx) = (1 - jTi 



2 /2 d"'P,(x) 



dx" 



Equations ([9|-([T0)i agree with the well-known formula for the associated Legendre functions 

(1-mV 
(l + m)\ 



(4) 
(5) 
(6) 



(7) 



(8) 



(9) 



(10) 



(11) 



(12) 



Different sign conventions appear in the literature, in particular regarding the place of the (- l) m , which is sometimes used in the 
definition of the associated Legendre functions instead (see, e.g., Chapter 8 of Abramowitz & Stegu n|1972| l. However, these sign 

O 

differences do not influence the vector spherical functions themselves. Also one sometimes considers l Ti m instead of T/ m as the 
toroidal vector spherical functions. 

From ( p~0| > and (|7]l-((8]l one immediately has 



T^„,(a,5) = (-l) m T* m (a,6), 
S l ,- m (a,5) = (-l) m S* lm (a,6). 



(13) 
(14) 
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2.2. Orthogonality properties 

The spherical functions Y\ m form an orthonormal sequence of functions on the surface of a sphere since 

f Y lm YJ, m , dQ = S"'5 mm ' , (15) 
Jn 

which is also complete in the Hilbert space S of the square-integrable functions on a sphere (making it an orthonormal basis of S 
in the L 2 norm). As in Fourier expansions, the completeness property is hard to establish and is related to the fact that spherical 
harmonics are eigenfunctions of a special kind of differential equations (see, for example, |Morse & Feshbach|1954[ ). This is accepted 
and not further discussed in this paper directed towards astronomical applications. Here and below one has dQ = cos 6 ddda, and 
the integration is taken over the surface of the unit sphere: < a < 2jt, -jt/2 < 6 < n/2, and <f J is the Kronecker symbol: S'-i - 1 
for i - j and <S y = 0, otherwise. 

Thanks to the completeness property a (square-integrable) complex-valued scalar function defined on a unit sphere f(a, 5) can 
be uniquely projected on the Yi m : 



oo / 



f(a,6) = Y J Y j fi m Y lm , 



1=0 m=-l 

where the Fourier coefficients /;,„ age given by 



flm 



= f fYl 
Jn 



dQ.. 



(16) 



(17) 



The equality in ( 16 1 strictly means that the right-hand side series converge (not necessarily pointwise at discontinuity points or 
singularities but at least with the I? norm) to some function /, and given the completeness and the orthonormal basis, the coefficients 
//,„ are related to / by ( 17 1. A truncated form of ( 16 1 with I < Z max < oo is an approximate expansion for which the equality in ( 16 1 
does not strictly hold (in general). 

Similarly, the VSH form a complete set of orthonormal vector functions on the surface of a sphere (with the inner product of the 
L 2 space): 



X 
X 
X 



W ctnm' 



T lm ■ r Vm , dQ = 6" 6' 



Sim ■ S* Vm , dQ = 6 U 6"' 



S [m -T; m ,dQ = 0. 

As we noted after (|3]l above, one also has the orthogonality between two vectors in the usual Euclidean space, 

Sim ' Tim = 



(18) 
(19) 
(20) 

(21) 



which holds for any point on the sphere. 

Any (square-integrable) complex-valued vector field V(a,6) defined on the surface of a sphere and orthogonal to u (radial 
direction), 



V(a, 6) = V a (a, 6) e a + V s (a, 5) e s , 

can be expanded in a unique linear combination of the VSH, 



/ 



V(a, 8) = } \ } t {timTim + SimSlJj , 



1=1 m=-l 

where the coefficients f/„, and s/,,, can again be computed by projecting the field on the base functions with 



Sin 



.. f v-t; 

Jn 

: r v-si 

Jn 



dQ, 



dQ. 



(22) 
(23) 

(24) 
(25) 
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2.3. Numerical computation of the vector spherical harmonics 

Analytical expressions of the VSH are useful for the lower orders to better understand what fields are represented with large-scale 
features and also to experiment by hand on simple fields in an analytical form. This helps develop insight into their properties 
and behaviour, which pays off at higher orders when one must rely only on numerical approaches. This is also a good way to test 
a computer implementation by comparing the numerical output to the expected results derived from the analytical expressions. 
Analytical expressions for the vector spherical harmonics of orders I < 4 are given explicitly as a function of a and 6 in Appendix 

m 

To go to higher degrees, one needs to resort to numerical methods. To numerically compute the two components of T/ m (a, 6) 
and Si, n (a, 6) given by |7])-([8]l, one first needs to have a procedure for the scalar spherical harmonics. Given the form of the Yi m (a, 6) 
in d9]l, the derivative with respect to a is trivial, and for S only, the derivative of the Legendre associated functions needs special 
care. There are several recurrence relations allowing the Legendre functions to be computed, but not all of them are stable for high 
degrees. Other difficulties appear around the poles when 5 + n/2 where care must be exercised to avoid singularities. Non-singular 
expressions and numerically stable algorithms for computing the VSH components are available in the literature, and the expressions 
we have implemented and tested are detailed in Appendix [5] 



2.4. Expansions of real functions 



Describing real-valued vector fields with complex-valued VSHs leads to unnecessary redundancy in the number of basis functions 
and fitting coefficients. In general, the coefficients in ( |23| l are complex even for a real vector field. But for a real field V(a, S), the 
expansion must be real. Given the conjugation properties of the VSH, it is clear that for a real-values function one gets 



ti m = (-V) m f l _ m , 
s lm = (-\) m s]_ m 



(26) 
(27) 



Then from the decomposition of the coefficients into their real and imaginary parts as 



s lm - s i„ t + 1 sf m ■ 



(28) 
(29) 



one gets at the end by rearranging the summation on m > 



F(a,<$) = ]T 



rji'K 

Im ' 



f 3 T 3 <R 
'1m 1 1m + S lm S lm J /m' 



S 3 ) 



(30) 



which is real. It is obvious that tf Q — sf = and, therefore, f/o = and sio = s^. 

It is useful to note that the use of real and imaginary parts of the VSHs as in ( 30 1 is analogous to the use of "sin" and "cos" 
spherical functions suggested e.g. by Brosche ( 1966| l. However, we prefer to retain notations here that are directly related to the 
complex-values VSHs and the corresponding expansion coefficients. 



Equation ( 30 1 provides the basic model (with only real numbers) to compute the coefficients by a least squares fitting of a field 
given for a finite number of points, which are not necessarily regularly distributed. In this case the discretised form of the integrals 



in ( 18 1, ( 19 1, and ( 20 1 is never exactly or 1, and the system of VSH on this set of points is not fully orthonormal. Then one cannot 



compute the coefficients by a direct projection. However, one can solve the linear model (30 1 on a finite set of coefficients f/ m and 
S[ m . Provided the errors are given by a Gaussian noise, the solution should produce unbiased estimates of the true coefficients. This 
is discussed further in Section 



2.5. Relation between the expansion in vector and scalar spherical harmonics 



Once we have fitted a vector field on the model ( 30 1, we have the expansion in the form ( 23 i with / < l max . The components V a and 
V s of the vector field V in the local basis e a and e$ as in ( 22 1 are expressed with the same set of coefficients f/,„ and si m on their 
respective components of the VSH. But since the V" and Vcue scalar functions on a sphere, they could also have been expanded 
independently in terms of scalar spherical functions Y\ m as V" = Tlm=-i ^im an d V 6 - Z™o Tlm=-i ^f m ^' m ' providing another 
representation, which has been used as mentioned in the Introduction to analyse stellar catalogues in the same spirit as with the 
VSHs. It is therefore important to relate the two different sets of coefficients for the purpose of comparing and discussing the major 
difference between the two approaches (see Section[3]). This issue is not central to this paper so the mathematical details are deferred 
to in AppendixO One must, however, note that the relations between coefficients Vf m and Vf , on the one hand, and t\ m and Si m , on 
the other, are ratner complicated and involve infinite linear combinations. It means, for example, that th e information contained in a 



single coefficient t/ m is spread over infinite number of coefficients V? and Vf (see, e.g. Vityazev 2010 for a particular case I = 1) 
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(a)5 u (b)5 2 ,i 



















fflpMj 









(c) S s 



(d) 5,0,5 



Fig. 2. Examples of spheroidal harmonics vector fields: S 14, $10,5, where one sees the change in the angular resolution 

with higher degrees. The S u is a glide flow in the direction of the y-axis. 



3. Transformation under rotation 

3.1. Overview 

A very important mathematical and practical feature of the expansion of a vector field on the basis of the VSHs is how they transform 
under a rotation of the reference axes. In short, when a vector field is given in a frame S and decomposed over the VSHs in this 
frame, one gets the set of components f/,„, s; m in this frame. Rotating the reference frame S (e.g. transforming from equatorial 
coordinates to galactic coordinates) to the frame S, one gets the transformed vector field that can be fitted again on a set of VSHs 
defined in the frame S. This results in a new set of components f;„, and si m , representing the vector field in this second frame. 
Significant additional work may be required to carry out the whole transformation of the initial field and perform the fit anew in 
the rotated frame. Fortunately due to the narrow relationship between VSHs and group representation, there is in fact a dramatic 
shortcut to this heavy procedure that adds considerably to the interest of using the VSH expansions in astronomy. This is illustrated 
in the self-explanatory accompanying commutative diagram, showing that to a space rotation (in the usual three-dimensional space) 
R corresponds to an operator % acting in the space of the VSHs and allows transforming t\ m and s\ m given in the initial frame into 
an equivalent set in the rotated frame: 

S — ^ S 

VSH 



{Urn, Sl,„} > \tlm,Slm} 

The amount of computation is negligible since this is a linear transformation between the coefficients of a given degree /. Overall, 
the set of t\ m and si m is transformed into itself without mixing coefficients of different degrees. 

This operator K has been introduced by E. Wigner in 1927 as the D-matrices. A convenient reference regarding its definition and 



properties is the book of Varshalovich, Moskalev & Khersonski ( 1988 1. Section 7.3 of that book is specifically devoted to this topic. 



The Wigner matrices have been used for decades in quantum mechanics, but probably they are not so well known in fundamental 
astronomy. We confine ourselves here to a quick introduction of the main formulas and conventions used in this paper. All the 
mathematical formulas have been implemented in computer programs (independently in FORTRAN and Mathematica), which can 
be requested from the authors. 

For practical astronomical applications, there is an unexpected benefit in this transformation under space rotation, which is not 
shared by the separate analysis on spherical harmonics of the components like Aa cos 6 and AS, or their equivalent for the proper 
motions. In general the two components for a given source obtained from observations do not have the same accuracy and are 
correlated. In the case of space astrometry missions like Hipparcos or Gaia, the correlations are generally smaller in the ecliptic 
frame due to the symmetry of the scanning in this frame compared to the correlations in the equatorial frame, where the solution 
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(a) T L 






(c) r 5 . 



(d) r 10 , 3 



Fig. 3. Examples of toroidal harmonics vector fields: 7\i, T^i, T$ ,3, Tio.s. One sees the change in the angular resolution with higher 
degrees. The 7\i is a simple rotation about the y-axis. 

is naturally available. Performing a least squares adjustment to compute f/„, and s\ m with correlated observations is a complication 
and standard pieces of software often allow only for diagonal weight matrix. One can easily get round this difficulty by carrying out 
the fitting in a frame where the correlations can be neglected (assuming such a frame does exist), and then rotate the coefficients 
with the Wigner matrix into the equatorial frame. This is equivalent to rigorously performing the fit in the equatorial coordinates 
with correlated observations. On the other hand, by fitting the components with scalar spherical harmonics, it is cumbersome to take 
the correlations between the components properly into account, and the coefficients in one frame do not transform simply into their 
equivalent into a rotated frame. 

More generally, if we allow for the correlations in the frame where the least squares expansion over the VSHs is carried out, 
one can still apply the Wigner rotation to the results and obtain the expansion in a second frame. The results would be exactly the 
same as if a new least squares fit had been performed in this frame by rigorously propagating the non-diagonal covariance matrix 
of the observations to this second frame to generate the new weight matrix. Using the Wigner matrix in this context is conceptually 
much simpler and in keeping with the underlying group properties of the VSHs, even though there might be no definite advantage 
for numerical efficiency (but no disadvantage as well). 



3.2. Mathematical details 

Consider two rectangular Cartesian coordinate systems S and S of the same handedness. Two coordinate systems are related by a 
rotation parametrised by three Euler angles a, b, and c. For a given vector with components x in S, the components x in S read as 



x = R z (c)R y (b)R z (a)x, 

where R, and R v are the usual matrix of (passive) rotation, respectively, about axes z and y 



( cose sine 0' 
■ sin e cos e 
1 J 



(31) 



(32) 



R, (e) ■ 



' cos e - sin e ^ 

1 
v sin e cos e 



where the order of arguments in the rotational matrices in ( 3 1 1 must be exactly followed. The inverse transformation reads 
x = R z (-a)R Y (-b)R z (-c)x. 



(33) 



(34) 
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Here the convention z - y - z is used for the sequence of rotations instead of the usual z — x-z that is more common in dynamics 
and celestial mechanics. The former choice has the great advantage of giving a real matrix d' nm (b) (see below) for the intermediate 
rotation, and this convention is universally used in quantum mechanics and group representation. The values of a, b, and c for 
transformations between the three usual frames in astronomy are listed in Table. [T] The relation between the z — y — Z and z — x — z 
conventions is given by 



R z (c) Ry(b) R z (a) = R z (c - n/2) R x (b) R z (a + n/2) , 
where R, are the usual matrix of rotation about axis x: 



Rx(c) : 



(10 \ 
cos e sin s 
- sin e cos s 



(35) 



(36) 



Then the transformations between the scalar and vector spherical functions defined in S and S are similar to each other and read 



as 



Y lm (a, 8) = D' m , m (a, b, c) Y lm ,(a, 5) , 

m'=-l 
I 

fi m {a, 5) = J] D l m ,Ja, b, c) T lm , (a, 5) , 

m'=-l 
I 

SiJa, 8)=^ D l m , m (a, b, c) S M {a, 5) , 



(37) 
(38) 
(39) 



rnf=—l 



where a and 6 are angles defined as right ascension and declination (generally speaking, the longitude and latitude of the spherical 
coordinates) from the components x, while a and 5 are those derived from the components x. In both cases Eq. Q is used. Here 
D l mn (a,b,c), \m\ < I, \n\ < I are the Wigner D-matrices (generalised spherical functions) defined as 



DUa,b,c)=e 



-i(ma+nc) jl 



(40) 



d l mn (b) = (-l) m ~" ^(l + m)\(l-m)\(l + n)\(l-n)\ 



"■niiix 



(cos I) 



2l—2k—m+n / 



Sill 



I) 



2k+m— n 



k=k a 



k\(l -m- k)\(l + n- k)\{m -n + k)\ 



&min = max(0, n - m), k max = min(/ - m,l + n) . 



(41) 



Many efficient ways to evaluate D'(a,b,c ) and d' mn {b) can be found in relevant textbooks. For example, Section 4.21 of 



Varshalovich, Moskalev & Khersonski 



( 1988 ) gives the expression of d l mn (b) in the form of a simple Fourier polynomial 



(42) 



where d' mn can be precomputed (using the obvious simplification of (41 ) as constants and used in further calculations for a 
given rotational angle b. However, even a straight implementation of (|40|-(4T i does the job without any problem for low to medium 
degrees. 

The inverse transformations of the spherical functions read as 



Y lm {a,8)= j] D l * m ,{a,b,c)Y lm ,{a,8), 

m'=-l 
I 

T !m (a, <5) = Z D mm'( a > b > c ) T im'(a, 5) , 

m'=-/ 
I 

Si m {a,8)= J] D l * m ,{a,b,c)S lml {a,8), 



(43) 
(44) 
(45) 



where the superscript '*' denotes complex conjugation. 

Now what matters for our purpose is the transformation of the VSH coefficients f; m and s/ m as in ( |23] l under a rotation of the 
coordinate system, and, more generally, that of coefficients fi m in the case of an expansion of a scalar field in the scalar spherical 
harmonics as given by ( |16) . 
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Given the completeness of the scalar and vector basis functions, we have defined the following expansions for scalar field / and 
vector field V on sphere (cf. ( 16 1 and (23 1): 



CO / 



CO / 



/ - / ; / l flm Y/m - Yj flm ^ lm ' 
1=0 m=-l 1=0 m=-l 

CO / CO / 

y = Yi ( f /m ^6n + S; m ) = ^ (f/„, Tfa, + 5/ m 5/,„) . 



(46) 
(47) 



1=1 m =-l 



1=1 m=-l 



Substituting (f3"7j)-( 39 1 or (|43]l-([45]l into (|46|-( 47 i and changing the order of summations, one gets the important transformation 
rules for the coefficients 



f lm = Yj Dl mm>( a , b , c ).flm> > 
m'=—l 
I 

tim= Y D l mm ,(a,b,c)ti m r , 



=-/ 



Sim* ■ 



S lm = Yj D 'mm'( a > b > C ) 

m'=-l 

and the inverse transformations 



/ 



fbn= Y D m'm( a < b ' C )flm' ■ 
m'=-l 
I 

tim= J] D ^ m (a,b,c)ti ml , 

m'=-l 
I 

Y D l :, m (a,b,c)s M . 



(48) 
(49) 
(50) 

(51) 
(52) 
(53) 



m'=-l 



Table 1. Euler angles a, b and c applicable to the change of astronomical frame. The obliquity has been taken for J2000 and 
equatorial frame can be seen as the same as ICRF at this level of accuracy. 



Transformation 


a 


b 


c 




deg 


deg 


deg 


Equatorial to Ecliptic 


270.0 


23.4393 


90.0 


Equatorial to Galactic 


192.8595 


62.8718 


57.0681 


Ecliptic to Galactic 


180.0232 


60.1886 


83.6160 



Once the expansion has been computed in one frame, it is very easy to compute the expansion in any other frame related by a 
rotation from the initial frame. Numerical checks have been performed with the analysis of a vector field in two frames, followed 
by the transformations (direct and inverse) of the coefficients f/ m and .?/,„) with (48]l-(|50]l and (|5T]l-(|53 I. 

These transformation rules allow us also to prove the invariance of the Eua 
degree: 



'Pi - Yj f lm flm - Yj f' m f lm ' 
m=—l m=—l 
I I 

7*1 = Y t,m tf'x ~ XI ~ tlm ~ t,m ' 

m=—l m=—l 
I I 

P\ - Yj S,m S lm ~ Yj * /m ' 



idean norm of the set of coefficients of a given 

(54) 
(55) 
(56) 



=-l m=-l 

This can be easily demonstrated by substituting (|48|-([50| into corresponding equation and using the well-known relation 



i 



YD l mm ,(a,b,c)D l : m „(a,b,c) = 6, 

n=—l 



Jm'm" ■ 



(57) 
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which is a special case of the general addition theorem for the Wigner D-matrices ( Varshalovich, Moskalev & Khersonski|1 988 
Section 4.7). The invariance of Pi will be used in Section 5.2 to assess the level of significance of the coefficients. 



4. Relation of the first-order VSH to global features 

The multipole (VSH) representation of a vector field has some similarity with a Fourier or a wavelet decomposition since as the 
degree / increases, one sees smaller details in the structure of the field. In the case of the comparison between two catalogues, 
increasing the degree reveals the zonal errors of the catalogue, while the very first degrees, say / = 1 and / = 2, show features 
with the longer "wavelengths". In particular, the first degree both in the toroidal and spheroidal harmonics is linked to very global 
features, such as rotation between the two catalogues or a systematic dipole displacement, like a flow from a source to a sink located 
at the two poles of an axis. It must be stressed that these global features are found in a blind way, without searching for them. They 
come out naturally in the decomposition as the features having the least angular resolution, and they are interpreted as a rotation or 
a dipole glide a posteriori. 



4.1. Rotation 



Consider a infinitesimal rotation on a sphere given by an the rotation vector R whose components in a particular frame are 
(Ri,Ri,R3)- This rotation generates the vector field 



V R = R x u , 

where u is the unit radial vector given by Q. Equation (58 1 can be projected on the local tangent plane and leads to 



V R = V R ■ e n 



-Ri sin 6 cos a - R2 sin 6 sin a + R3 cos 6 , 



V R = V R -e s = R l sina-R 2 cosa. 



(58) 



(59) 
(60) 



From the explicit expressions of the toroidal harmonics given in Table A. 1 one sees that any infinitesimal rotation field has non-zero 



projection only on the T\ m harmonics, and is orthogonal to any other set of toroidal or spheroidal harmonics. Therefore (58 1 can be 
written as a linear combination of the T\ m . One has explicitly the equivalence 



V R - ffo^io + ffi^n + h.-i^i-i 



where 



j? = _jt* 
'11 




1 R 2 ), 



(61) 

(62) 
(63) 



allowing the three parameters R of the rotation to be extracted from the harmonic expansion of degree I — 1 . A change of reference 
frame will lead to the same rotation vector expressed in the new reference frame. Given the importance of a rotation between 
two catalogues or in the error distribution of a full sky astrometric catalogue, having the rotational vector field as one of the base 
functions is a very appealing feature of the VSH. This contrasts with the analysis of each components of the vector on scalar spherical 
harmonics, in which a plain rotation projects on harmonics of any degree as shown in | Vityazev] ( | 1 997| [20 1 ) or in Appendix [C] 
However, the analysis of each component as scalar fields retains its interest if one has good reasons to assume that the field has been 
generated by a process in which the two spherical coordinates have been handled more or less separately. In this case there is even 
no reason to change the frame for another one since the two components may have very different statistical properties. 



4.2. Glide 

The rotation field is the most common fully global signature on a spherical vector field, and often it even seems to be the only one 
deserving the qualification of global. However, if we consider a rotational field G xu, one may draw at each point on the sphere a 
vector perpendicular to the rotation field and lying in the tangent plane to the sphere. In this way we generate another field V c with 
axial symmetry, fully orthogonal to the rotation field G xu. This new field has no projection on the T\ m harmonics and, in fact, on 
no other T/ m . Its components are given by 



V£ = -Gisina + G 2 cosar, (64) 
Vg = —G\ sin 6 cos a - G 2 sin 6 sin a + G3 cos 5 , (65) 

where the vector G = {G\,G 2 , G3) determines the orientation and the magnitude of the field, as R did for the rotational field. Going 



to Table A. 1 it is easy to see that it is equivalent to the S\ m field with the decomposition 



V 6 = i? Sio + iftSn + if^Si-i (66) 
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where 

(67) 

4 ^(-G 1+ iG 2 ), (68) 

showing its close relationship with a rotation field. A corresponding picture can be seen in Fig.|2]for an axis along the y-axis. 

The same can be stated somewhat differently: we have introduced the glid^jvector field V, which is orthogonal to both radial 
vector u and the rotation field G x u, or equivalently one has 

V c = hx(Gxh) = G-(Gh)h, (69) 

which shows clearly that the associated vector field is just the projection of G on the surface of the unit sphere (we have above G 
minus its radial component), that is to say the projection of an axial, or dipole, field on the sphere. 

As a pattern, a glide field is as global a feature as a rotation, and can be seen as a regular flow between a source and a sink 
diametrically opposed. From the astronomical point of view this is a field associated to a motion of the observer toward an apex, 
with all the stars showing a kinematical stream in the opposite direction. For extragalactic sources like quasars, this also characterises 
the effect induced on QSO systematic proper motions resulting from the acceleration of the observers with respect to the frame where 
the QSOs are at rest on average fFanselow|1983| |Sovers et a/. 1 1 998] |Kovalevsky |2003 [ Kopeikin & Makarov 2006 ; Titov, Lambert 
|& Go ntier 2011). The main source of this acceleration is thought to be the centripetal acceleration from the galactic rotation. Its 
unambiguous determination should be achieved with Gaia. It is important to recognise that in the general analysis of catalogues 
or spherical vector fields, rotation is not more fundamental a feature than glide, and both must be searched for simultaneously. 
Otherwise, for most of the distributions of the data on the sphere, when the discrete orthogonality conditions between the VSHs 
are not fully satisfied, there will be a projection of the element not included in the model in the subspace generated by the other 
one. Thus, fitting only the rotation, the rotation vector will be biased if there is a glide not accounted for in the data model. And 
reciprocally for the glide without rotation properly fitted. 



(a) Rotation 
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(b) Glide 



Fig. 4. Rotation and glide with the z-axis shown here with a Mercator projection to emphasise the close relationship between the 
two vector fields. The constant displacement in longitude from the rotation has its counterpart with a constant glide displacement in 
the reduced latitude (Mercator latitude). 



4.3. Duality of the rotation and glide effects 

Finally the duality between the infinitesimal rotation and glide can be made striking if the fields are plotted in the Mercator pro- 
jection, as illustrated in Fig. |4] Both effects are along the z-axis and have the same magnitude for R and G. Thus they are rotation 
and glide with the same axis. The effect of a small rotation is seen as a uniform translation in longitude, identical to all latitudes, 

1 We coined this word to have something as easy to use as rotation to name the transformation and the feature and to convey in a word the idea 
of a smooth flow between the two poles. 
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as expected for a cylindrical projection. For the glide, the translation is uniform in the reduced latitude given by the Mercator 
transformation 

6 M = In (tan + £)) (70) 

and does not depend on the longitude. 

There is a clear similarity between the two fields (both represent a mapping of the sphere onto itself with two fixed points: 
the poles), but also an important difference seen with the Mercator projection: for the rotations, the mapping of a flow line (the 
continuous line tangent to the vectors of the field) is bound between the two extreme boundaries -n and +n, which must be identified 
as the same line on the sphere. The rotation axis is the natural source of a cylindrical symmetry shown by the Mercator projection 
(or any other cylindrical projection). For the glide, it goes very differently since the two end points of a flow line are sent to infinity, 
in both directions, and can only be identified in the projective space, a much more complex manifold than the Euclidean cylinder. 
Could this important topological difference be the underlying reason why only the rotation is usually considered as the only obvious 
global structure? We do not have the answer and leave this issue open for the moment. As long as infinitesimal transformations are 
concerned, which is exactly what matters for the analysis of astrometric catalogues, the two effects are in fact very similar and must 
be taken into account together as a whole. 



5. Practical implementation 

In principle, once the difference between two catalogues has been formed or when a vector field on a sphere (e.g. a set of proper 
motions) is available, it is normally trivial to project the underlined vector field on the VSH, by computing the integrals (24i 
with an appropriate numerical methods. This will work as long as the distribution of the sources on the unit sphere preserves the 
orthogonality conditions between the different base functions. What is precisely meant is the following: using the same numerical 



method to evaluate the Fourier integrals, one should check that Eqs. ( 1 8 I, ( 19 1 and ( 20 1 hold with a level of accuracy compatible with 
the noise up to some degree / max . Unless one has many sources almost evenly distributed, this condition is rarely fulfilled. However, 
the fact that the coefficients can be expressed by integrals, that is to say by direct projection on the base functions, is a simple 
consequence of the application of a least squares criterion leading to a diagonal normal matrix when the orthogonality condition 
holds. Starting one step back to the underlying least squares principle we can directly fit the coefficients to the observed components 



of the field with the model ( 30 1, by applying a weighting scheme based on the noise in the data. In case the orthogonality conditions 



hold, the solution will be equivalent to the evaluation of the integrals. When this assumption is not true, meaning with the discrete 
sampling and/or uneven weight distribution at the source coordinates, the harmonics are not numerically orthogonal, but the least 
squares fitting provides an unbiased estimate of the coefficients t\ m and si m , together with their covariance matrix. 



5.1. Least squares solution 



The fitting model is given by ( 30 1 up to the degree / = Z max . It is fitted to a set of N data points sampled at (a k , 6 k ), k — 1 , 2, • ■ • , N, 



where at each point the two components of the vector field A (e.g., difference between two catalogues, or proper motion components) 
are given as 



A(a k ,6 k ) 



A a k e a + A 6 k e 6 , 



1,2, 



,N, 



(71) 



together with the estimates of their associated random noise o"a« and cr A «. For the sake of simplicity we assume here that the 

k k 

covariance matrix between the observations is diagonal, meaning the observations are realisations of independent random variables. 
The observations are fitted to the model ([30| with the least squares criterion 



mm 

ttm.Slm 



S 

k=l 



' (A« - V a (a k , 6 k )f (A* - V s {a k , S k j) 



2\ 



A? 



(72) 



(*», t 



From (30 1 and (22 1 we see that with N observations one has 2N conditional equations to determine the 2/ max (/ max + 2) coefficients 



sf m ),m=l,2,...,U=l,2,...,Z n 



l /m' Vm' Sl0 ' s lm> s im)' "i = 1, 2, ...,/,/= 1, 2, ... , / max . The size of the design matrix is 2N x 2/ max (Z max + 2), which can be very 
large in analysis of astrometric surveys, with N above 10 5 - 10 6 and a sensitivity on a small enough angular scale to represent the 
regional errors leading to explore up to / max ~ 15 — 20. The data storage could become a problem if one wishes to store and access 
the full design matrix. One solution to get round this problem is to build up the normal matrix on the fly, taking its symmetry into 
account. For a large data set (N » 'm ax X this bookkeeping part will be the most demanding in terms of computation time. For fixed 
N the computational time for this part of calculations grows like f mzx . In general the problem is rather well conditioned, with small 
off-diagonal terms, and is not prone to generate numerical problems to get the least squares solution. 



5.2. Significance level 

The least squares solution ends up with an estimate for each of the coefficients included in the model. The model must go to a certain 
maximum value of Z max by including all the component harmonics of order m, for each I < Z max . The total number of unknowns is 
then 2Z max (/ max + 2) and grows quadratically with Z max . Within a degree /, one can decide on the significance of each coefficient of 
order m, from its amplitude compared to the standard deviation. However, these individual amplitudes are not invariant by a change 
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of reference frame, like the components of a vector in ordinary space, and the significance of a particular component is not preserved 
through a rotation of the frame. On the other hand, the whole set of f/,„ or s/ m for a given / behaves like a vector under rotation, and 
the significance (or lack of) of the vector itself, that is to say its modulus, is intrinsic and remains true or false in any other rotated 
frame. Therefore it is more interesting and more useful in practice to investigate the significance level for each degree I rather than 
that of the individual coefficients. The reason behind this statement is deeply rooted in the invariance property of the Y[ m under a 
rotation in the three-dimensional Euclidean space and the fact that the set of spherical harmonics F/ m of degree I is the basis of an 
irreducible representation of the rotation group of the usual three-dimensional space. Therefore, this set is globally invariant under 
a rotation, meaning there is linear transformation linking the Y\ m in one frame to their equivalent of same degree / in the other frame 
(see Section[3]l. 

Therefore, one should consider the 21+1 elements of the set F/ m , for m = -/,.../ like the components of a vector that are 
transformed into each other under a rotation, keeping the Euclidean norm unchanged in the process. Given the linearity of the 
derivative operator, this property extends naturally to the T\ m and S/ m . For the same reason, the coefficients f/,„ and S\ m in the 
expansion of a vector field must be considered together within a degree I and not examined separately, unless one has good reasons 
to interpret physically one or several coefficients in a particular frame: the Oort coefficients in the Galactic frame are one of the 



possible examples, as shown in Vityazev ( 2010 1. More generally a wide range of parameters in Galactic kinematics have very natural 
representations on vector harmonics, such as the solar motion towards the apex or a more general description of the differential 
rotation than the simple Oort's modelling (see for example Mignard (2000|). 

Therefore it is natural to look at the power of the vector field and at how this power is spread over the different degrees I. Let P 
be the power of a real vector field V(a, 5) on a sphere defined by the surface integral)^] 

P(V)= f \V\ 2 d£l, (73) 
Jn 

or its discrete form used hereafter with numerical data, 

nv)^ — Ywf, (74) 

n 

i=\ 

where n is the number of samples on the sphere. 



When projected on the VSH one finds easily with ( 18 I, ( 19 1, and (23 1 that 



i=i 

m=-Z 
/ 

*?=j>6»£.- (75) 



For a real vector field having the symmetries ( 26 1— ( 27 1 one gets 

i i 



I I 

n = 4 + 2 ^ \s lm \ 2 = 4 + 2 £ ( (sf m f + (sfj) . (76) 



m=l m=l 



The quantities P\ and PJ are the powers of toroidal and spheroidal components of the degree I of the field. Both quantities are scalars 
invariant under rotation of the coordinate system (see Section[3jl. Thus P' t and P s { represent an intrinsic property of the field, similar 
to the magnitude of a vector in an ordinary Euclidean vector space. 

We can now take up the important issue of the significance level of the subset of the coefficients s\ m and f/ m for each degree /. In 
harmonic or spectral analysis, this is always a difficult subject, because it needs to scale the power against some typical value that 
would be produced by a pure noisy signal. In general one can construct a well-defined test, with secure theoretical foundations, but 
with assumptions about the time or space sampling, which are rarely found in the real world. Therefore the theory provides a good 
guideline, but not necessarily a fully safe solution applicable in all situations. The key issue here is to construct a criterion that is 
well understood from the statistical view point and robust enough to be usable for mild departure from the underlying assumptions, 
like regular time sampling in time series. We offer below three possibilities for testing whether the power P\ and PJ is significant 
and therefore that there is a true signature of the VSH of degree I in the vector field. In the following we drop superscripts V and 
'f' and write simply Pj. Relevant formulas below are valid separately for P' t and PJ or for P\ + P*. The involved coefficients will be 
denoted as r;,„ to represent either f/„, or si„, or both. 



2 The power is often defined with a normalisation factor 1/47T before the integral. The present choice leads to simpler expressions for the power 
in terms of the coefficients f;„, and s/ m and has been preferred. 
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Start from the null hypothesis Hq that the signal is just made of noise, producing some values of Pi in the analysis. One needs 
to build a test to tell us that above a certain threshold, the value of Pi is significant and would only be produced by pure noise 
with a probability p <k 1. To build the test properly, the probability density function (PDF) of Pi is required. It can be obtained 
with the additional assumptions that the sources are rather evenly distributed on the sphere. The least squares fitting also gives us 
the variance-covariance matrix of the unknowns si, n and t\ m . Given the orthogonality properties of the VSH, the normal matrix is 
almost diagonal if the distribution of sources over the sky and the associated errors of the quantity representing the vector field are 
homogeneous. In this case the covariance matrix is also almost diagonal. 

Consider the coefficients r®, rf\ and r? in ( 30 1 with a fixed degree I. Each diagonal term of the relevant part of the normal 
matrix (unweighted here) takes the form 



for 
for 
for 



no, 

'lm ■■ 



(77) 
(78) 
(79) 



where the sum extends over all the sources on the sphere. Here R[ m stands for either 7/ m or S /,„. The first term with m = does not 
depend on the right ascension a, while the others have a part in cos 2 mor and sin 2 ma, whose average is 1 /2 for a ranging from to 2n 
and a relatively regular distribution in a. Given the normalisation of the VSHs, we now see that the diagonal terms with m + are 
just two times the term with m = 0. The covariance matrix being the inverse of the normal matrix, one eventually has the variances 
of the unknowns as 



var(r, ) = cr 



/o ■ 



var(O=^ /2, 

' lm' 



var(r^)=cr? /2, 



where crf Q is the variance of r/Q. Then dividing Pi in (76 1 by cri. leads to the normalised power: 



W, 



Pi 

^2 



z 



+ (^ 



(80) 
(81) 
(82) 



(83) 



With ( 8 1 1-( 82 1 and the assumption that the signal is pure white Gaussian noise, this is the sum of 21 + 1 squares of standard normal 
random variables with zero mean and unit variance. It is well known that W/ follows a x 2 distribution with n — 21 + 1 degrees of 
freedom. With this normalisation, W/ is strictly proportional to the power Pi. 

Under the Hq hypothesis that the signal is just white noise, one can build a one-sided test to decide whether a harmonic of degree 
/ is significant with the probability level y: 



P(W, > w) 



X 2 2l+1 (x)dx = 



r(/+ 1/2, w/2) 

r(/ + 1 12) 



(84) 



where F(a, x) and T(a) are the incomplete and complete gamma functions, respectively. This quantity can be easily computed. In 
some cases it may be easier to use the transformation of W ilson & Hilferty| ( |1931 1, 



Z = 



(?f-K) 



(85) 



where n is the number of degrees of freedom. It is well known that Z approximately follows, even for small n, a standard normal 
distribution of zero mean and unit variance. A test risk of y — 0.01 is achieved if Z > 2.33 or at the level of y — 0.025 is Z > 1.96. 
Experiments on simulated data have shown that the key elements in the above derivation, namely that cr /0 =s cr/ m ,m > is 
practically achieved when sources are rather regularly distributed in longitudes, even with irregularities in latitudes. The reason is 
basically that this factor V2 follows from the average of cos 2 ma and sin 2 ma, while the latitude effects are shared in a similar way 
in all the r; m . 



Alternatively, instead of scaling the power by only cr 2 () and taking advantage of the asymptotic properties ( 80 1-( 82 1, it would 



have been more natural to define a reduced power with 



^te)+z 



W" lm> 



cr 



2\ 



/;;/ I 



(86) 



which is not strictly proportional to the power. But what matters in practice is whether one can construct an indicator, related to the 
power, with well defined statistical properties, so one may relax the constraint of having exactly the power, provided that a reliable 
test of significance can be built. For least squares solutions with Gaussian noise, and small correlations between the estimates of 
the unknowns, it is known that each coefficient^follows a normal distribution of zero mean whose variances can be estimated with 
the inverse of the normal matrix. Therefore, Wi is also a sum of squares of standard normal variables and its PDF is that of a x 1 
distribution with 21+1 degrees of freedom. For regular distributions of the data points, the two reduced powers W/ and W/ are 
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equivalent. So either test should give the same result. In ( |86) the multiplying factor two before the sum in f76] > has been absorbed 
by var( \2Y) = 2 var(F) (i.e. by erf- and erf being approximately a factor V2 smaller than <x ;o ). 



Finally, we can also go back to the first test in (83 I, which is the power scaled by the variance of the component siq. Since 



the variance aj is in practice an estimate of a random variable based on the data, it may happen that the value produced by 
the least squares is unusually small, which will make W/ too large and trigger a false detection. Therefore, while this test has a 
good statistical foundation for regular distribution of sources, it is not very robust to departure from this assumption. We suggest 
improving this robustness, a very desirable feature in the real data processing, by replacing the scaling factor cri. by the average of 

the set ( <t 2 , 2(cr^) 2 , 2(crf m ) 2 ), m = I,.,., I whose mathematical expectation is the true value of cr 2 , but it has a smaller scatter than 



cr 2 Let cr 2 be this average, we have then instead of ( 83 i a new expression for the reduced power 



Wi = Wi = -f , (87) 

which is directly proportional to the power and follows ax 2 distribution with 21+1 degrees of freedom. 

From a theoretical standpoint, the three indicators are equivalent for a large number of data points evenly distributed on a sphere. 
In practical situations, with some irregularities in the angular distribution and/or a limited number of sources, the significance testing 
with Wi should by construction be more robust against random scatter affecting the estimated standard deviations. Preliminary 
experiments have confirmed this feature, although a wider Monte Carlo simulation is needed to quantify the advantage. 



6. Application to FK5 Catalogue 

The comparison between the FK5 and Hipparcos was published soon after the release of the Hipparcos catalogue by Mignard & 
Froeschle (2000) primarily provides the orientation and the rotation (spin) of the FK5 system with respect to the Hipparcos frame. 
The differences in position and proper motions are shown in Fig.[5]for the whole set of stars used in the comparison. In the positional 
plot one clearly sees a general shift toward the south celestial pole and few regional effects of smaller scales. 

The global analysis carried out at that time was restricted to fitting these six parameters from the positional and proper motion 
differences from Hipparcos to enable the transformation of astronomical data from the FK5 frame into the Hipparcos frame or, 
more or less equivalently, to the ICRF. The remaining residuals were large and considered as zonal errors, not reducible to a simple 
representation and presented in the form of plots of Aacos 6 and A6 as a function of the position of the stars. Large zonal errors 
were clearly shown by these plots, but no attempt was made to model these residuals. 

Using the VSH now, with the software developed during this work, it is possible to draw refined conclusions and to show that the 
terms of low degree beyond the rotation are truly significant. In particular we can see in Table [4]that the positional glide term from 
the spheroidal harmonic of degree I = 1 is larger than the orientation term. Terms of degree / = 2 are also large, both in position and 
proper motion, a feature that had not been noticed as clearly with the earlier analysis based on a priori model. 




(a) Positions (b) Proper motions 



Fig. 5. Differences in positions and proper motions between the FK5 catalogue and Hipparcos in 1991.25. The typical arrow sizes 
in the plots are respectively 100 mas and 2.5 mas/yr. 



6. 1 . Global differences 

This concerns the rotation and the glide between the two systems, or equivalently using the VSH terminology, the toroidal and 
spheroidal harmonics with I - \. The results for the rotation shown in Table [2] are fully compatible with the earlier analysis of 
Mignard & Froeschle (2000) and confirm a systematic difference between the FK5 inertial frame and the ICRS materialised by 



Hipparcos. The interpretation of the systematic spin is also discussed in the same paper and seems to be related to the precession 
constant. These new values also agree with a similar analysis performed by |Schwan| (|200r), provided the epoch is corrected from 
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1949.4 (mean epoch of the FK5 observations) to 1991.25 (Hipparcos Catalogue reference epoch). A new result is given in Table[3] 
for the glide, but is only significant in the positions. It describes a systematic difference of about 60 cos 6 mas in the declination 
system between the two frames. This is a very large systematic effect given the claimed accuracy of the FK5 below 50 mas. Nothing 
similar is seen for the proper motions. Any transformation of a catalogue tied by its construction to the FK5 system and later on 
aligned to the HCRF (Hipparcos frame aligned to the ICRF) must absolutely take this large systematic difference into account in 
addition to the rotation. The results just shown correspond to the harmonics I = 1 from a fit of the position and proper motion 
difference to l m . dx = 10. They are slightly different with another choice of the maximum degree, but within the formal error. In both 
tables the standard deviations quoted are derived from the covariance matrix of the least squares solution scaled by the r.m.s. of the 
post-fit residuals. 



Table 2. Global orientation (in 1991.25) and spin between the FK5 and the Hipparcos catalogue. 
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0.09 


e ; : 


18.5 


2.4 


o) z : 


0.82 


0.09 



Table 3. Glide between the FK5 and the Hipparcos catalogue in positions (in 1991.25) and proper motions. 





Position (mas) 
J1991.25 


a 




PM(mas/yr) 

cr 


gx- 


18.3 


2.4 


7x- 


0.04 0.09 


Sy- 


-1.3 


2.4 




0.18 0.09 


Sz- 


-64.0 


2.4 


7z- 


-0.37 0.09 



6.2. Regional differences 

The analysis at higher degree, up to / = 10, clearly shows significant power in most degrees, indicating a definite structure in the 
space distribution of the differences. Given the absence of zonal error in the Hipparcos catalogue, at the level seen in this analysis, 
they must come entirely from the FK5. A detailed analysis of this regions variations was performed by Schwan (2001) by projecting 
separately each components (difference in right ascension and declination) on a set of scalar orthogonal polynomials, instead of using 
the vectorial nature of the differences. Given the numerous positional instruments used to produce the FK5, some measuring only 
the right ascension, others only the declination, others giving both, analysing the components in equatorial coordinates separately 
makes sense and may reveal structure scattered in many degrees with the VSHs. Here we simply want to illustrate the relevance 
of the VSH to analyse the difference between two real catalogues. A detailed discussion of the application to the FK5 system is 
deferred to an independent paper. 



6.2.1. Regional differences in position 



The amplitudes, defined as the square root of the unweighted power, in each degree for the positional differences between the 
FK5 and Hipparcos in 1991.25 are given in Table [4] for the toroidal (left part of the table) and spheroidal (right part of the table) 
harmonics. The level of significance computed with Eq.(87 1 and the normal approximation of Eq. (85 i is shown in the columns 
labelled Z. With I = 1 one recovers the very large significance of the rotation and the glide. Harmonics I = 2 are also very large, 
explaining together with the glide a significant part of the systematic differences between the two catalogues. Then the powers 
remain statistically very significant up to / = 7, although with small amplitudes, and then decline. With the values of the // m and s\ m 



(not given here), one could produce an analytical transformation relating the two systems, in the same spirit as in Schwan (2001 



With the last column one sees that the expansion to / = 10 explains about 58% of the variance seen in the data for the position and 
32% for the proper motion. 



6.2.2. Regional differences in proper motions 

The same analysis has been done for the proper motions. The powers and their significance level expressed with a standard normal 
variable are listed in Table [5] Again the signal is strong in the harmonics of the first two degrees and up to / = 7 primarily for the 
toroidal components. 
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Table 4. Amplitude (power 1 / 2 ) and significance level based on (87) for the positional differences between FK5 and Hipparcos 
expanded on the VSH. 



degree 


I 2 \ 1/2 
\Xm tj,,,) 


Z 


1 2 \ 1/2 
\2jm S l m ) 


Z 


Rem. (power) 1 ' 2 




mas 




mas 




mas 













474.5 


1 


85.8 


10.18 


192.8 


19.90 


425.0 


2 


96.6 


11.48 


139.5 


15.94 


389.6 


3 


53.6 


6.01 


52.8 


5.89 


382.3 


4 


65.6 


7.44 


66.2 


7.52 


370.8 


5 


96.8 


11.41 


62.2 


6.72 


352.5 


6 


65.2 


6.95 


82.2 


9.36 


336.5 


7 


68.3 


7.18 


42.2 


3.00 


326.7 


8 


45.0 


3.18 


57.1 


5.22 


318.6 


9 


31.1 


0.28 


41.4 


2.25 


314.3 


10 


34.4 


0.63 


50.6 


3.63 


308.3 



Notes. The last column gives the power in the remaining signal after the VSH model up to degree / 
gives the square root of the power of the data as defined by l|73[>-l|74|l. 



has been removed. The first line with 1 = 



Table 5. Amplitude (power 1 / 2 ) and significance level based on (87) for the proper motion differences between FK5 and Hipparcos 
expanded on the VSH. 



degree 


(Zm hm) 

mas/yr 


Z 


(Zm S/,,,) 

mas/yr 


Z 


Rem. (power) 1/2 
mas/yr 













14.6 


1 


3.1 


9.85 


1.2 


3.67 


14.2 


2 


2.3 


7.40 


2.5 


8.00 


13.8 


3 


1.1 


2.10 


1.5 


3.88 


13.7 


4 


1.9 


5.46 


1.9 


5.47 


13.4 


5 


2.7 


7.97 


1.3 


2.34 


13.1 


6 


2.0 


5.16 


1.8 


4.33 


12.8 


7 


2.8 


8.02 


1.3 


1.44 


12.4 


8 


1.2 


0.52 


1.5 


2.26 


12.3 


9 


1.3 


1.17 


1.2 


0.28 


12.2 


10 


1.7 


2.91 


1.2 


-0.04 


12.0 



Notes. The last column gives the power in the remaining signal after the VSH model up to degree Z has been removed. The first line with / 
gives the square root of the power of the data as defined by l|73[>-l|74[l. 



7. Analysis of the expected Gaia results 

Another important application of the VSH technique is an analysis of mathematical properties (e.g. systematic errors) of an astro- 
metric catalogue, as well as extraction of physical information from the catalogue. Here we concentrate on the QSO catalogue as 
part of the expected Gaia catalogue. Our goal is to investigate (i) the anticipated accuracy of the determination of the rotational state 
of the reference frame and the acceleration of the solar system with respect to the QSOs, and (ii) the expected estimate of the energy 
flux of the primordial (ultra-low-frequency) gravity waves. From the mathematical point of view, this amounts to checking the accu- 
racy of determining the VSH coefficients of orders / = 1 and / = 2. The toroidal coefficients of order 1 describe the rotational state 
of the reference frame with respect to the QSOs, while the spheroidal coefficients of order 1 give the acceleration of the solar system 
with respect to the QSOs, which shows up as a glide ( |Fanselow|198 3 ; Sover s~ef q/.|1998| Kovalevsky |2003 1 |Kopeikin & Makarov| 
2006; [Tito v, Lambert & Gontier|201l) . The VSH coefficients of order 2 can be related to the energy flux of the ultra-low-frequency 
gravity waves ( Gwinn et al.\1997\ . 



7.1. Simulated QSO catalogue from Gaia 

We simulated the Gaia QSO catalogue with all the details that we know at the time of writing. The total number of QSOs observable 
by Gaia - 700000 - was chosen according to the analysis of Mignard ( 2012|l. The realisti c distribution of the QSOs in the Gaia 
integrated magnitude G was taken from Slezak & Mignard ( 2007)> (see also |Robin et at 2012|l. The distribution of the Gaia-observed 
QSOs over the sky was chosen according to the standard Gaia extinction model ( [Robin et a/.|2012| l. The model follows the realistic 
distribution of the dust in the Galaxy, so that the QSO distribution used in this section is not homogeneous on the sky and depends 
on both angular coordinates with the most pronounced feature being the avoidance area close to the Galactic plane. The proper 
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motions of QSOs are generated from a randomly selected rotation of the reference frame and an acceleration of the solar system's 
barycentre with respect to the QSOs with magnitude and direction expected from the circular motion of the solar system with respect 
to the galactic centre (i.e. the proper motions are generated from the full set of VSH coefficients of order 1). These resulted in our 
simulated "true" (noise-free) catalogue of QSOs. 

The expected accuracies of astrometric parameters (positions, proper motions, and parallaxes) as functions of the G magnitude 
and the position on the sky are taken from the standard Gaia science performance model (ESA 2011 de Bruijne 2012). Then 
the QSO simulated solution with Gaia is obtained by adding Gaussian noise with the corresponding standard deviation to each 
parameter of each source. 

One more feature used in our simulations is that we take the correlations between the estimates of the two components of proper 
motion into account. As a model for the correlations, the empirical distribution is calculated from the corresponding histogram for 
the Hipparcos catalogue as shown on Fig. 3.2.66 of ESA ( 1997 vol. 1, Section 3.2). We neglect herewith the complicated distribution 
of the correlation on the sky as shown on Figs. 3.2.60 and 3.2.61 of ESA|( [T997] l. The generated random correlation is then used to 
generate correlated Gaussian noise for the components of the proper motion of each source. 



7.2. Main results of the analysis of the simulated QSO catalogue 

From the resulting realistic QSO catalogue expected from Gaia, we fitted the VSH coefficients for maximal order Z max ranging 
from 1 to 15. The correlations between the components of the proper motion for the same star are taken into account in the fit by 
using block diagonal weight matrix with 2x2 blocks corresponding to each source. Accounting for these correlations in the fit 
generally modifies the estimates at the level of one tenth of the corresponding standard deviations. We repeated these simulations 
(including generation of the QSO catalogue) tens of times. For a typical simulation, the resulting biases of the estimates of the 
rotation of the reference frame and the acceleration of the solar system, as well as the standard errors of those estimates are shown 
in Fig. [6] The biases are computed as differences between the estimated values and true ones (those put in the simulated 'noise-free' 
catalogue). For the standard errors of the estimates, we took the formal standard deviations from the fit multiplied by the square root 
of the reduced^ 2 of the fit. Comparing of the estimates for different values of / max as shown in Fig. [6] is a useful indicator of the 
reliability of the estimate. Our analysis generally confirms the assessments of the anticipated Gaia accuracy for the reference frame 
and acceleration published previously: both the rotational state of the reference frame and the acceleration of the solar system could 
be assessed with an accuracy of about 0.2 yuas/yr. However, this accuracy can only be reached under the assumption that (1) Gaia 
will be successfully calibrated down to that level of accuracy (that is, the errors of astrometric parameters remain Gaussian down 
to about 0.2 yt/as), and (2) real proper motions of the QSOs (transverse motions of the photocentres due to changing QSO structure, 
etc.) are not too large and do not influence the results. Whether these assumptions are correct will not be known before Gaia flies. 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

I™* 



(a) Rotation 



1 2 3 4 5 6 7 



9 10 11 12 13 14 15 



(b) Acceleration (divided by the light velocity) 



Fig. 6. True errors and formal errors of the estimates of the absolute values of the rotation and acceleration from the fits with 
1 < /max ^ 15 for the simulated QSO catalogue matching expected Gaia accuracy. 



Finally, we discuss the expected upper estimate of the energy flux of the ultra-low-frequency gravity waves that Gaia could 
measure. We use the computational recipe described by |Gwinn et al. ( 1997| see Eq. (11) and the discussion in the first paragraph 
of their Section 4). Our analysis shows that Gaia can be expected to give an estimate at the level of Oqw < 0.00008 IT 1 for the 
integrated flux with frequencies v < 6.4 x 10 9 Hz. Here, h = i//(100km/s/ Mpc) is the norma lised Hubble constant. This is to 

( 1997 i. The same approach applied to 
1.5 x 10~ 9 Hz. Therefore, one can 



be compared with the VLBI estimate Qr.w < 0. 1 1 /z 2 for v < 2 x 10 9 Hz of |Gwinn et al. 



the VLBI catalogue used by Titov, Lambert & Gontier (2011 



_ gives Q GW < 0.009 fr 2 for v < 
expect that Gaia will improve current estimates from VLBI data by two orders of magnitude while covering a larger interval of 
frequencies. However, this conclusion relies on the ability to see large-scale features of small amplitude in the Gaia proper motion 
thanks to the improvement of the detection in l/NJ, UTCes . This can be hindered by correlated noise in the intermediates results or 
more likely the various sorts of unmodelled systematic errors, not known today, although every effort is made in the instrument 
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design and manufacturing to keep them low and compatible with the above results. These systematic errors can be expected to 
influence lower degree VSH coefficients to a greater extent than the higher degree ones. Therefore, for analysing the gravity waves, 
it is advantageous to use not only the quadrupole harmonics as in Gwinn et al. ( \991) and|Titov, Lambert & Gontier (2011), but 
also an ensemble of harmonics that includes higher order harmonics (see e.g. |Book &"F lanagan 201 l| l. Given the number of QSOs 
in the Gaia catalogue, such a refined approach is feasible and will help establish more reliable limits on the energy flux resulting 
from primordial gravitational waves. 



7.3. Correlations between the VHS coefficients 

It is interesting to briefly discuss the correlations between the VSH harmonics resulting from the least squares fits with the simulated 
catalogue of 700 000 QSOs described above. These correlations were computed form the inverse of the weighted normal matrix from 
the VHS fit. Depending on the maximal order / max , the root mean square value of the correlations is 0.04, while a few correlations 
may attain the level of 0.4. Having such a large number of sources one might expect that the correlations should be extremely 
small. Indeed, for a catalogue of 700000 sources homogeneously distributed on the sky, the r.m.s. of the correlations between the 
VSH coefficients is typically 0.0006 and the maximum is at the level of 0.006. We see that the correlations for the realistic QSO 
catalogue are about 70 times larger. One can find two reasons for the larger correlations: (1) the inhomogeneous distribution of 
the QSOs on the sky (this can be called the "kinematical inhomogeneity" of the catalogue) and (2) the fact that the accuracies of 
the QSOs' proper motions, which are used in the weight matrix, also depend on the position on the sky (this can be called the 
"dynamical inhomogeneity" of the catalogue). For the QSO catalogue used in this section the dynamical inhomogeneity alone gives 
correlations with an r.m.s. of 0.01 and the maximum of 0.2, while the kinematical inhomogeneity alone leads to the r.m.s. of 0.03 
and the maximum of 0.35. Generally speaking, correlations lead to higher errors in the fitted parameters. We finally note that the 
dynamical inhomogeneity will be slightly greater for the real Gaia catalogue than we see in our simulation. This is caused by the 
fact that the standard Gaia science performance model (ESA 2011 |de Bruijne|2012 i is smoothed in the galactic longitude. However, 



our calculations show that this additional inhomogeneity does not significantly change the results given above. 



8. Conclusion 

In this paper we have shown the relevance of the vector spherical harmonics (VSH) to decomposing, more or less blindly, spherical 
vector fields frequently encountered in astronomical context. This happens very naturally in comparing of stellar positional or 
proper motion catalogues or the different solutions produced during the data analysis in global astrometry. We have provided the 
mathematical background and considered several practical issues related to the actual computation of the VSH. 
We have in particular: 

- shown how to perform the decomposition with a least squares minimisation, enabling processing of irregularly distributed sets 
of points on the sphere; 

- provided a way to test the significance level of any degree in the expansion and then set a criterion for where to stop the 
expansion; 

- stressed the importance of the invariance properties against rotation of the expansion and given the necessary formulas that 
express the decompositions in the most common astronomical frames; 

- provided applications with a reanalysis of the FK5 catalogue compared to Hipparcos, and an example of how this tool should be 
used during the Gaia data processing to prepare the construction of the inertial frame and derive important physical parameters; 

- compiled an extensive set of practical expressions and properties in the Appendices. 
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Appendix A: Explicit formulas for the vector spherical functions with I < 4 



Table A.l. Toroidal harmonics T [m (a, 6) and spheroidal harmonics S; m (a, 5) for 1 < I < 4. 

Harm. Mult. Components 

coef. e a e 6 

T m hJk cos 6 



T 



1 1 



S21 
S22 



S 



S 



2 \ 2jr 

I ^| sin 5 (cos a + 1 sin a) - sin a + l cos a 



r 20 i^g sin2o 

i - cos 2<5 (cos a + 1 sin a) - sin 6 (sin a - 1 cos a) 

T 2 2 I - sin 25 (cos 2a + l sin 2a) 2 cos 5 (sin 2a - l cos 2a) 

r 30 \^ cos 5 (5 sin 2 6 - 1) 

r 3] ^ ^ sinc5(15 sin 2 <5 - 11) (cos a + I sin a) -(5 sin 2 5 - 1) (sin a - i cos a) 

j\ - cos 6 (3 sin 2 6 - 1 ) (cos 2a + 1 sin 2a) sin 26 (sin 2a - 1 cos 2a) 

^33 ~k \ ~ co&2 $ sm $ ( cos 3a + 1 sin 3a) - cos 2 o (sin 3a - 1 cos 3a) 







r « ]| J| sin 2(5 (7 sin 2 (S - 3) 

^ 41 fa (28 sin 4 o - 27 sin 2 o + 3) (cos a + 1 sin a) - sin 6 (7 sin 2 (5 - 3) (sin a - 1 cos a) 

7,42 il \l - sin 2<5 (7 sin 2 (5 - 4) (cos 2a + 1 sin 2a) cos 6 (7 sin 2 6 - 1) (sin 2a- 1 cos 2a) 

^ 43 16 yl cos 2 <5 (4 sin 2 <5 - l)(cos3a+ 1 sin 3a) -3 cos 2 S s'm6 (sin 3a- 1 cos 3a) 

^ 44 \yh - cos 3 sine) (cos 4a + 1 sin 4a) cos'o (sin 4a- 1 cos 4a) 



2 y 2jr 

\ yjl sin a - 1 cos a sin 6 (cos a + 1 sin a) 



S w \J± cos<5 



4 \ 2jr 



S20 7 sin2o 



I ^| sin 6 (sin a - 1 cos a) - cos 26 (cos a + 1 sin a) 

I yj^ -2 cos o (sin 2a - 1 cos 2a) - sin 26 (cos 2a + 1 sin 2a) 

\^ coso(5sin 2 o- 1) 

S31 th sf^ (5sin 2 o- l)(sina- 1 cosa) sino (15 sin 2 <5 - 1 1) (cos a + 1 sina) 

^ 32 I y^c - sin 2o (sin 2a - 1 cos 2a) - cos 6 (3 sin 2 6 - 1) (cos 2a + 1 sin 2a) 

^33 Tf, v ~ cos 2 o (sin 3a - 1 cos 3a) cos 2 <5 sin 6 (cos 3a + 1 sin 3a) 



±Jj sin2<5(7sin 2 <5-3) 

^ 41 ~k -\fi sind(7sin 2 o- 3) (sin a- 1 cosa) (28 sin 4 o - 27 sin 2 o + 3) (cos a + 1 sina) 

S « Tesfl ~ cos 6 ( 7 sin2(5 ~ 1 ) ( sin 2a ~ * cos 2a ) ~ sin 2(5 ( 7 sin ^ ~ 4 ) ( cos 2a + * sin 2ar ) 

^ 43 TE sjl 3 cos 2 (5 sin (sin 3a - 1 cos 3a) cos 2 (5(4 s'm 2 6 - 1) (cos 3a + 1 sin 3a) 

^ 44 I y^L - cos'o (sin 4a - 1 cos 4a) - cos 3 <5 sin <5 (cos 4a + 1 sin 4a) 
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Appendix B: Practical numerical algorithm for the scalar and vector spherical harmonics 

B.1. Scalar spherical harmonics 

Scalar spherical functions F; m can be computed directly using definitions (|9|>— (|T0]> . The Legendre functions can be evaluated numer- 
ically with the following stable recurrence relation on 



(/ - m) Pim(x) = (21 - V)xPi-i M (x) -(l-l+m) P ; _ 2 , ; 
starting with 

Pfn— 1 , m — , 

P mm = (2m-\)\\(\-x 2 T 12 . 



(B.l) 

(B.2) 
(B.3) 



The algorithm based on (B.l i-( |B3] i is discussed in Section 6.8 of Press et til. (19921. The algorithm described there aims at 
computing a single value of Pi, n (x) for given values of I > 0, < m < I and \x\ < 1 . It is trivial to generalise it to compute and store 
all the values of Pi m (x) for / < Z max , Z max > 0, and a given x. The algorithm takes the form 



^00 = 1, 

Pm+l,m+l = ( 2m + 1) Vl - X 2 P mm ,m = 0, . . 
Pm+l,m ~ (2ffl + 1) X Pmm ? — 0, . . . , /max 

1 



l, 



Phn = 



I — m 



((21 - V)xPi-i, m -(1-1+ m)Pi- 2 „) , 



l = m + 2, 



,/n 



0, 



,/n 



■2. 



(B.4) 
(B.5) 
(B.6) 

(B.7) 



Here Z max > is the maximal value of I to be used in the computation. The notation "a = . . . , b — . . ." denotes outer cycle for b and 
inner cycle for a with the specified boundaries. As a result one gets a table of values of Pi m (x) for a specified x (\x\ < 1) and for all 
/ and m such that < m < I and < I < / max . An extension of this and related algorithms useful for very high orders / is given by 
|Fukushima| ( |20T2l ). 



B.2. Vector spherical harmonics 

To compute vector spherical functions as defined by — <[8j one also needs to compute derivatives of the associated Legendre 
functions. A closed-form expression for the derivatives reads as 



dPlm(x) 

dx 



Ix I "4" /// 

Plm(x) + ~ ^Pl-l tm (x). 



l-X 



1 



(B.8) 



Special care should be taken for the computations with 6 close to ±n/2. Indeed, for 6 = ±n/2, one has x = sin 6 = ±1, and both 
the factors 1/ cosi5 in <|7]l-([8]) and the derivative P' n (x) go to infinity. Therefore, numerical computations (in particular, Eq. (B.8i) 
become unstable for 6 close to +n/2. To avoid this degeneracy the definitions of T\ m and 5; m can be rewritten as 



Ti m (a, S) = (-!)" 



S lm (a, 5) = (-!)" 



21+1 (l-m)\ 



^ An 1(1+ 1) (l + m)\ 

x\Ai m (sinS)e a - I S te (sin <5) e s j 



21+1 (I -m)l 



\ An 1(1+1) (l + m)\ 

x ^i B fa (sin 6) e a + A; m (sin 5) 
A lm (x)= Vl -x 2 P' (x), 



B bn (x)=m 



vr 



: Plm(x). 



(B.9) 



(B.10) 
(B.ll) 
(B.12) 



One can easily see that (e.g. from (Hi) that Ai„,(x) and Bi„,(x) remain regular for any < m < I and \x\ < 1. It is easy to see that 
numerically stable algorithm for A/,,, and B/,,, read as 

5/0 = 0, / = l,...,/max, (B.13) 

B u = l, (B.14) 
_ (2m + l)(m + 1) r— 1 

L>m+l,m+\ — VI X tS mm , 

m 

w=l,...,Z max -l, (B.15) 
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\,m — (2/71 + 1 ) X B mm , fit — 1 , . . . , / ma x 1 ■ 

l 



B, 



I — m 



((21 - 1) x fl M -(/- 1 + m) B,^ m ) , 



l = m + 2,..., Z max , m = 1, . . . , Z max - 2 . 

A/ = Vl -X 2 B n , 1= l,...,/max, 

Afa = - (-i I B/,„ + (/ + m) B/_i, m ) , 
m 

m = 1,...,/, / = 1,.. . J m ax ■ 



(B.16) 

(B.17) 
(B.18) 

(B.19) 



As a result tables of A;„,(jc) and B/ m (x) are computed for a given |x| < 1, < m < I and 1 < / < / max . Explicit expressions for the 
vector spherical harmonics for 1 < I < 4 are given in Annex [A] 



Appendix C: VSH expansion of a vector field vs. the scalar expansions of its components 

As explained earlier with the VSH, we can expand a vector field defined on a sphere on vectorial basis functions preserving the 
vectorial nature of the field and behaving very nicely under space rotations. However, based on the usual practice in geodesy and 
gravitation, following Brosche ( 1966 1 and Brosche ( 1970| l astronomers have also used, for many years, the expansions of each 
component of the vector field, namely V a — V ■ e a and V = V ■ eg in the scalar spherical harmonics. It is interesting to show the 
relationship between the two expansions and how their respective coefficients are related to each other. 



C.1. Definition of the expansions 

Let us consider a vector field V on the surface of a sphere. Its VSH expansion reads as 

oo / 



V(a, 8) = ^ ^ (t lm T lm + s tm SiJj . 



(CI) 



/=1 m=-l 



On the other hand, the components of this vector field V, V s are also expandable independently of each other in terms of the usual 
scalar spherical functions Y/ m : 



CO / 



/=() m=-l 
oo I 

;=o m=-l 
or for the vector field itself: 

oo / 

n«.*) = Z Yj{ V ?>n e « + V ?>« e s) Yl " 

1=0 m=-l 



(C.2) 
(C3) 



(C4) 



As with all our infinite sums, the equalities (C.l i-(C.4i only hold in the sense of convergence in the L 2 norm as already stated in 
Section [2~2] Pointwise convergence is not guaranteed in the general case. 
Therefore, one has the identity 



2 X { t,mTlm + SlmSlm ) = Z Z K»' e ° + V l'm' e *) Y ''> 



(C5) 



1=1 m=-l 



l'=0 m'=-l' 



C.2. Formal relations between the coefficients 



Multiplying both sides of (C.5 1 by T* lm or S* lm (again, superscript '*' denotes complex conjugation) and integrating over the surface 
of the sphere, one gets 



CO /' 



tlm ~ ho ^\ ^ ( V " m ' A ' ml ' m ' ~ V '' m ' Blml ' m ) ' 



CO /' 



S[m ~ ^1(1 11 X i V " m ' Blml ' m ' + V '' m ' A ' ml ' m ) ' 

y V / /'=() m'=-V 



(C6) 
(C7) 
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where 



r dY\ m 

Mml'm' = I -T^T Yy m i dQ, (C.8) 



r 1 oy: 

B, m ,> m >= -—^Y hn ,dQ, (C.9) 

J n cos 5 da 

where as usual dQ. = cos 5d5da, and the integration is taken over the surface of the unit sphere: < a < 2n, -n/2 < 6 < nj2. 
On the other hand, multiplying both sides of (C.5 i by Y*, and integrating over the surface of the sphere one gets 

VS, = 2 J] ^ (t Vm ,A* Vmllm + SVm ,-B> Vm ,^, (CIO) 

= 1 m'=-l> V 1 V + V > 

OO I' ^ 

Vf m = J] £ tfnr+vi {- tVm ' B m*» + Sl ' m ' A i'^) ' (C11) 

l'=\ m'=-l' V v / 

It remains to compute A// m </ m and B/' m '/ m in a convenient way. But formally we have obtained the two-way correspondence between 
the coefficients (f /m , s /m ) and (V^, Vfj. 

C. 3. Explicit formulas for A hnVm , and Bi mVml 
It is straightforward to show that 

1 



Hml'm' 



'brd'i 



2 nd O yim7l+2k+l,m 

x (-la, mk + + m)B,_ UmMl ) , (C.12) 
— ml 7t 5 mm ' 6 l+2kJ ' y lm y l+2k , m B lmk , (C. 1 3) 



(l-m)l 

7l m =J(2l+l) y - J - (C.14) 

V (/ + m) ! 



-f 

* J-i 



X 

ai mk =- I f /+2*+i,m(*) fhW dx , (C.15) 

J-i Vl -x 2 

1 f 1 1 

Plmk=- , Pl+2k,m( x ) Plm(x) dx , (C.16) 

* J-i Vl - X 2 

where 6'i is the Kronecker symbol = 1 for i = j and d'i = otherwise) and k e Z is arbitrary integer. Both o-/,,^ and Bi mk are 
positive rational numbers provided that the indices m, and k are selected in such a way that the Legendre polynomials under the 
integral are different from zero: |m| < / and \m\ < I + 2k + 1 for ai mk , and |m| < / and \m\ < I + 2£ for /?/,„/;. Otherwise a\ mk and are 
zero. Numbers a\ mk and can be computed directly or by using recurrence formulas that can be obtained from the well-known 
recurrence formulas for the associated Legendre polynomials. The following relations allow one to consider only the case of B\ mk 
with positive indices: 

(Xljn-k = a l-2k+\,m,k-\ , (C.17) 
Pl,m,-k = Pl-2k,m,k > (C. 1 8) 

(l-m)\ (l-m + 2k+l)\ 
(1 + m)\ (I + m + 2k+ 1)1 
(Z-m)! (l-m + 2k)\ 

Pl,-m,k - 77- — rr 77- ; ^rrr Plmk , (C.20) 

(/ + my. (I + m + 2ky. 

aimk = 2/ + 4^ + 3 ( W ~ m + 2k + Z)fiijnJc+r 

+ {l + m + 2k+l)B lmk ). (C.21) 

Finally, a number of equivalent formulas for B\ mk valid for any / > 0, < m < I and k > can be derived. Thus, using that any 
associated Legendre polynomial can be represented as a finite hypergeometric polynomial and integrating term by term, one gets 
the following formula for Bi mk valid for any I > 0, < m < I and k > 0: 

_ (2m -1)!! 2 «-™ +k) (2m + 2p-W 
P,mk ~ 4f* 4o 2P(p + 2m)\ 
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mm(p,l-m) 



Z 



■m+s r,m+p-s 



l+2k 



s=m r dx(0,p-l+m-2k) 



sl(p-s)\ 



(a + b)l 
(a-b)lbl 



(C.22) 
(C.23) 



We note that a number of explicit formulas for j3i, n k can be found for /, m, and k satisfying some specific conditions (e.g. m = I). 
For the general case, an alternative formula for can be derived using the Gaunt formula (Ga unt|1929| Appendix, pp. 192-196). 
However, Eq. (C.22i is sufficient for practical computations. These formulas can be easily implemented numerically. Equation 
( |C.22| i involves a double sum of terms of alternating signs. This means that computational instabilities can appear if floating point 
arithmetic is used. 



C.4. Final relations between the coefficients 



Since A/ m /< m * and B/,„/< m < vanish unless certain relations between the indices hold, one can significantly simplify the transformations 
( C.6 l-( |C.7| i, and ( |C. 10] >— ( (C. 1 1 1 can be significantly simplified and represented as a single sum over integer k: 

< _ n W>» V (v a 

/_j \ v l+2k+ 1 ,m Plmk 



2y/l(l+l) 



k><- 



+ V\ 



l+2k,i 



Jt+2k,mPlmk) , 

1 /TTTTTT Zj \- V l + 2k,m ml n+2k,m/3lmk 

2 V'(' + 1) TT._, 



(C.24) 



^7 In 



k>'- 



+ Vf+2k+l,m Pi™*) , 



(C.25) 



v: 



lm 



2 



Z 



1 



/<>- 



V/ + 2k + 1 



tl+2k+l,m <]lmk 



^7 In 

2 



z 



m i yi+2k,m 

+ s l+2k,m , Plmk 

yfl + 2k 
1 ( mi yi +2 k. 



(C.26) 



/<>- 



V/ + 2k + 1 



-tl+2k 



VTh2^ 



/mA- 



Plmk - yi+2k+\,m (~l »lmk + (I + m) fil-\,m,k+\ ) 
7;+2/t+l,m 



?6nt - 



yjl + 2k + 2 



(-(/ + 2Jt+l) 



+ (/ + 2^ + 1 + m)j3i mk ) . 



(C.27) 
(C.28) 

(C.29) 



In the above formulas we simplified the notations of the lowest values of k in the sums and assume now that V" { = 0, Vf j = 0, 
f()m = 0, and so™ = and that the terms containing them vanish identically. 

Thus, we have proved that the information of each particular f/ m and si m is distributed over infinite number of V? and Vf, and 
vice versa. 



C.5. Relation to space rotations 



It is important here to draw attention to two major differences between the expansions ( C. 1 1 and ( C.2 1-( C.3 1. 



- A rotation between two catalogues is very simply represented in the vectorial expansion with the harmonic T\ m and the three 
coefficients ?i o, fi -i, fi,i, while from (C.lOi- (C.lll, this will require an infinite number of coefficients in the component 
representations. One may rightly argue that the opposite is equally true: a simple scalar decomposition only over, say, I = 1 will 
be much more complex in the vectorial expansion. This is true, but physically the really meaningful global effects are precisely 
the rotation and the glide, which generate very simple vector fields and project on VSHs of first degree only and requires many 
degrees if done with the scalar components. We know of nothing equivalent to the scalar representation, unless, obviously, one 
builds an ad-hoc field for the purpose of illustration. 
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- The behaviour of the VSHs or of the scalar spherical harmonics Y\ m under space rotation are very similar and show the same 
global invariance within a given degree I. Their transformations are fully defined with the Wigner matrix as shown in Section[3] 
Therefore it seems at first g l ance that the transformations under space rotation of the ti m and s/,„ in ( C. 1 1 are not simpler than that 
of the Vf and Vf m in (C.2 i— ( C.3 1. There is, however, a very important difference that makes the use of the VSH so valuable. By 

~~ , the new coefficients correspond to the expansion of the vector 



applying the Wigner matrix associated to a space rotation to ( C. 1 



field with its components given in the rotated frame. Now in (C.10i-(C.ll i, the two scalar fields are considered independently 



of each other, and the new coefficients deduced from the application of the Wigner operator correspond to the expansion of the 
initial scalar fields (V a , V s ) expressed in the rotated fields, but these components are not those of the field V projected on the 
rotated coordinates, and they are still the components in the initial frame, since a scalar field is invariant by rotation. 
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